Kronig-Penney Eigenenergies

Below is the expression for one part of the expression for the  Kronig-Penney eigenergies:

temp1=Sum[(1/(ϵ-((((k+2π n)))^2))),{n,-Infinity,Infinity}]//FullSimplify

Here I have copied that expression and made it into a function.   I then plot it for a given choice of k.  Note that the sharply sloped line passing through ϵ=40 is spurious - it's only there because Mathematica is trying to connect different diverging sides of the cotangent function.

f[k_,ϵ_]:= ((-Cot[(1/2) (k-(Sqrt[ϵ]))]+Cot[(1/2) (k+(Sqrt[ϵ]))])/(4 (Sqrt[ϵ])))

Plot[f[.2,ϵ],{ϵ,0,100},PlotRange-> {-1,1}]

To get the eigenengergies, we need to plot a horizontal line equal to 1/v0, and look for all points where it crosses the curve.  We can get a set of solutions for differnt k by plotting contours of the function f.

ContourPlot[{f[k,ϵ]==.5,f[k,ϵ]==2,f[k,ϵ]==1,f[k,ϵ]==0},{k,-π,π},{ϵ,0,50},PlotPoints-> 100]

Just for fun, here's the function f as a 2D plot.   The bands are obtained by slicing it at z=1/v0.

Plot3D[temp1,{ϵ,0,100},{k,-π,π},PlotPoints-> 100,PlotRange-> {-.1,.1}]

Let' s  calculate the curvature of lowest band.   We need to get an expression for ϵ(k) for small k and small ϵ.

temp3=Normal[Series[temp1,{k,0,2},{ϵ,0,0}]]

soln=Solve[temp3== (1/v0),ϵ]

Grab the first root, and expand as a function of k.

temp4=ϵ/. soln[[1]]
temp5=PowerExpand[Normal[Series[temp4,{k,0,2}]]]

If you look closely, you'll see that this is proportial to (k^2). This means that  we have an effective mass, since the energy is parabolic in k.   We can plot the approximate effective mass as a function of v0.  We'll see that as the potential gets bigger, the effective mass gets bigger.

Plot[((((( (1440 v0+240 (v0^2)+4 (v0^3)))/(10 v0 (((12+v0))^2)))))^(-1)),{v0,0,100}]